unit uGlsTree;
(*
  Description:
  Implementation of the tree class that generates a tree based
  on the modelparameters described in a TModelParameters class.
  A description of the model used here can be found in:

  Creation and Rendering of Realistic Trees,
  Jason Weber and Joseph Penn,
  SIGGRAPH 95 Conference Proceedings, Annual Conference Series,
  pp. 119-128, Addison Wesley, August 1995

  A link to the paper:
  http://www.cs.duke.edu/education/courses/fall01/cps124/resources/p119-weber.pdf
  Original source code on Delphi by Ton van den Heuvel, 2002
*)

(*
  The trunk needs a base (Bole) that can split and have Flare.
  Stems need to be able to Clone and Spawn
  Leaves need a Stalk, Pattern of placement, Direction (hanging down, layered; facing sun, upturned for rain)
  Fruit need " "
*)

interface

uses
  System.SysUtils,
  Vcl.Forms,
  System.Classes,
  Vcl.Dialogs,

  uModelParameters,
  uGlsFrame,
  uGlsONB,
  uGlsMatrix,
  uGlsVector,
  uTransformationList;

type
  PLeaf = ^TLeaf;

  TLeaf = record
    RotationX: real;
    RotationY: real;
    RotationZ: real;
    Origin: TPoint3d; { TPoint3d is declared in GlsVector }
  end;

type
  PFruit = ^TFruit;

  TFruit = record
    RotationX: real;
    RotationY: real;
    RotationZ: real;
    Origin: TPoint3d;
  end;

var
  // 5 defines of shapes-Layers, Points, (x,y,z point)
  TreePolygonVertices: array of array of array of real;
  TreePolyImageVertices: array of array of array of Integer;
  TreePolygonVerticesNumber: array of Integer; { # of points in each define }

  // Only 1 Fruit shape.. no defaults
  FruitPolygonVertices: array of array [0 .. 2] of real;
  FruitPolyImageVertices: array of array [0 .. 1] of Integer;
  FruitPolygonVerticesNumber: Integer;

  // Leaves are wasting space: should be Oval and User Defined
  { const } { Leaf definitions, hard coded for the moment...later on these leaf
    definitions should be loaded from leaf definition files. }
  { oval, thin, triangle, round, 3-lobe maple,
    3-5-lobe maple, 5-lobe maple, (5-lobe oak,) 3 leaflets,
    User Defined  This is just a start... }
  { 9:0..8(types)  ?(points per leaf)  3:0..2  (x,y,z) }
  LeafPolygonVertices: array of array of array of real;
  LeafPolyImageVertices: array of array [0 .. 1] of Integer;
  { array[0..0, 0..8, 0..2] of real =
    (((0, 0, 0), (-0.2, 0.1, 0), (-0.3, 0.5, 0), (-0.2, 0.9, 0),
    (0, 1, 0), (0.2, 0.9, 0), (0.3, 0.5, 0), (0.2, 0.1, 0),
    (0, 0, 0))); }
  LeafPolygonVerticesNumber: array [0 .. 1] of Integer { = (9) };
  { for i := 0 to LeafPolygonVerticesNumber[LeafShape] do
    begin
    LocalPoint.x := LeafPolygonVertices[LeafShape, i, 0]; }

type
  TStem = class; // Predefined so it can be itself

  TTree = class
  private
    FTrunk: TStem; { Need a List or Array of Stems... Base Split }
    FModelParameters: TModelParameters;
    FTrunkLength: real;
    FScaleTree: real;
    FTrunkRadius: real;
    FTotalStems: Integer;
    { FTotalTreeLeavesGrown, }
    FTotalGrownLeaves: Integer;
    FTotalFruit: Integer;

    FMinX, FMaxX: real; // Variables holding the maximum x, y and z values
    FMinY, FMaxY: real; // occurring in the tree
    FMinZ, FMaxZ: real;
    FSeason: byte;

    // FQuality sizes and changes # of Leaves grown
    { ?? specifies the amount of points a section is made up of }
  public
    property TrunkLength: real read FTrunkLength;
    property ScaleTree: real read FScaleTree;
    property TrunkRadius: real read FTrunkRadius;
    property TotalStems: Integer read FTotalStems write FTotalStems;
    property TotalGrownLeaves: Integer read FTotalGrownLeaves
      write FTotalGrownLeaves;
    { property TotalTreeLeavesGrown: integer read FTotalTreeLeavesGrown write FTotalTreeLeavesGrown; }
    property TotalFruit: Integer read FTotalFruit write FTotalFruit;
    property MaxX: real read FMaxX write FMaxX;
    property MinX: real read FMinX write FMinX;
    property MaxY: real read FMaxY write FMaxY;
    property MinY: real read FMinY write FMinY;
    property MaxZ: real read FMaxZ write FMaxZ;
    property MinZ: real read FMinZ write FMinZ;
    property Trunk: TStem read FTrunk;
    property Season: byte read FSeason;
    property ModelParameters: TModelParameters read FModelParameters;

    constructor Create(const AModelParameters: TModelParameters;
      const ASeason: byte);
    destructor Destroy; override;

    procedure Grow;
    { Pre: FModelParameters <> nil
      Post: A tree using the model parameters specified in
      FModelParameters has been created. }
  end;

  TStem = class
  private
    FSections: TList;
    FBaseLength: real;
    FLength: real;
    FRadius: real;
    FTree: TTree;
    FSubStems: TList;

    FLeafs: TList;
    FFruit: TList;
    FParent: TStem;
    FLengthChildMax: real; // The maximum relative length of the substem
    FOffsetChild: real; // Holds the current y position along the stem
    // in global coordinates
    FSectionFrameList: TList; // List of frames of the sections that make up
    // the current stem. This list when we start
    // spawning child stems from the current stem
    FStemClones: TList;
    FLengthCloneMax: real;
    FOffsetClone: real;
    FCloneFrameList: TList;
  public
    property BaseLength: real read FBaseLength;
    property Length: real read FLength;
    property Radius: real read FRadius;
    property Tree: TTree read FTree;
    property Sections: TList read FSections;
    property SubStems: TList read FSubStems;
    property Leafs: TList read FLeafs;
    property Fruit: TList read FFruit;
    property LengthChildMax: real read FLengthChildMax;
    property OffsetChild: real read FOffsetChild;
    property SectionFrameList: TList read FSectionFrameList;
    property StemClones: TList read FStemClones;
    property LengthCloneMax: real read FLengthCloneMax;
    property OffsetClone: real read FOffsetClone;
    property CloneFrameList: TList read FCloneFrameList;
    constructor Create(const ATree: TTree);
    destructor Destroy; override;

    procedure CreateStructure(const AParent: TStem;
      const ALength, AOffsetChild: real; const Level: Integer);
    { Pre: True
      Post: The complete tree structure has been created. The total number of
      stems and leaves in the tree is known. The length of the stem
      is needed because the amount of substems that will be spawned from
      the current stem depends on the length of the current stem. }
    { createStem    createBranch    splitStem   spawnBranches }

    { ZDistribution getBranchDistribution
      ( int pass, double offset, parlen, length ); }
    function CalculateBranchDistribution(const Level: Integer;
      const y, StemLength: real): real;

    function CalculateSectionRadius(const Level: Integer;
      const y, StemLength, StemRadius: real): real;
    { Pre: 0 <= y <= 1
      Ret: The radius of the stem at the (normalized) y position along the stem.
      See for the exact calculations the paper "Creation and Rendering of
      Realistic Trees", by Jason Weber and Joseph Penn. }
    function CalculateVerticalAttraction(const Level: Integer;
      const SectionFrame: TCoordinateFrame): real;
    { Pre: ALevel > 1 (there is no vertical attraction for the trunk and main branches)
      Ret: An angle in radians that is added to the segments curve angle to simulate
      vertical attraction. }
    function ShapeRatio(const Shape: Integer; const Ratio: real): real;
    { Pre: 0 <= Shape <= 8 & 0 <= Ratio <= 1
      Ret: A certain predefined ratio depending on the Shape and Ratio parameter.
      For precise details see "The Creation and Rendering of Realistic Trees". }

    procedure Grow(const AFrame: TCoordinateFrame; const ARadius: real;
      const Level: Integer);
    { Pre: AFrame <> nil
      Post: A stem at the recursion level specified by ALevel has been created,
      with length ALength and maximum radius ARadius, using
      coordinate system AFrame, growing on tree TTree. }

    procedure GrowSubStems(const Level: Integer);
    { Pre: Level >= 0
      Post: The sub stems for the current stem have been created. }
    procedure GrowStemCellClones(const Level: Integer);
    { Pre: Level >= 0
      Post: The  stems Clones OF..for the current stem have been created. }
    procedure GrowLeaves(const Level: Integer);
    { Pre: Level >= 0
      Post: The leaves for the current stem have been created. }
    function AdjustLeafFrame(const ALeafFrame: TCoordinateFrame)
      : TCoordinateFrame;
    procedure GrowFruit(const Level: Integer);
    { Pre: Level >= 0
      Post: The Fruit for the current stem have been created. }

    procedure CreateSection(const AFrame: TCoordinateFrame; SectionRadius: real;
      const Vertices: Integer);
    { Pre: AFrame <> nil
      Post: A list of points that makes up the created section has been added
      to the FSections list. }

    procedure CreateLeaf(const LeafShape: Integer;
      const AFrame: TCoordinateFrame);
    { Pre: AFrame <> nil
      Post: A list of points that makes up the created leaf has been added to
      the FLeafs list. }
    procedure CreateFruit(const AFrame: TCoordinateFrame);
    { Pre: AFrame <> nil
      Post: A list of points that makes up the created Fruit has been added to
      the FFruit list. }
  end;

  TSection = TList;

implementation

uses Math,
  uConstants, {Globals}
  fGlsProgress, {POV Calls to build tree}
  fGlsTreedemo, {Calls to build tree}
  fGlsViewer {Calls to build tree};

function GetRandomValue(const UpperBound: real): real;
{ Pre: True
  Ret: A random value i, with i: -UpperBound <= i <= UpperBound }
var
  Sign: Integer;
begin
  if UpperBound = 0 then
    Result := 0
  else
  begin
    { Determine 'nondeterministically' the sign of the result }
    if Random(2) = 0 then
      Sign := -1
    else
      Sign := 1;
    { Calculate the result using Random that returns a real value i, i: 0 <= i < 1 }
    { Multiply by < 1 is dividing so result is -UpperBound.. UpperBound }
    Result := UpperBound * Random * Sign;
  end;
end;

{ TTree }
constructor TTree.Create(const AModelParameters: TModelParameters;
  const ASeason: byte);
begin
  inherited Create;
  { Clones creates the LEVEL.. not input any where else }
  FModelParameters := AModelParameters.Clones;
  FTrunk := TStem.Create(Self);
  FMaxX := 0;
  FMaxY := 0;
  FMaxZ := 0;
  FTotalStems := 0;
  { FTotalTreeLeavesGrown:=0; }
  FTotalGrownLeaves := 0;
  FTotalFruit := 0;
  FSeason := ASeason;
  with FModelParameters do
  begin
    FScaleTree := Scale + GetRandomValue(ScaleV);
    FTrunkLength := (nLength[0] + GetRandomValue(nLengthV[0])) * FScaleTree;
    FTrunkRadius := FTrunkLength * Ratio * Scale0 { * BaseSplits ? };
  end;
end;

destructor TTree.Destroy;
begin
  FTrunk.Destroy;
  FModelParameters.Destroy;
  inherited Destroy;
end;

procedure TTree.Grow;
{ Pre: True
  Post: A tree using the model parameters specified in FModelParameters has
  been created. }
var
  AFrame: TCoordinateFrame;
begin
  { Create a standard canonical frame }
  AFrame := TCoordinateFrame.Create;

  { Create trunk }
  with FModelParameters do
    FTrunk.Grow(AFrame, FTrunkRadius, 0);
end;

{ TStem }
constructor TStem.Create(const ATree: TTree);
begin
  inherited Create;
  FSections := TList.Create;
  FSubStems := TList.Create;
  FStemClones := TList.Create;
  FLeafs := TList.Create;
  FFruit := TList.Create;
  FSectionFrameList := TList.Create;
  FTree := ATree;
end;

destructor TStem.Destroy;
var
  i, j: Integer; // Loop variables
begin
  { Destroy all segments in the segment list }
  for i := 0 to FSections.Count - 1 do
    for j := 0 to TList(FSections[i]).Count - 1 do
      TPoint3d(TList(FSections[i])[j]).Destroy;

  for i := 0 to FLeafs.Count - 1 do
    for j := 0 to TList(FLeafs[i]).Count - 1 do
      TPoint3d(TList(FLeafs[i])[j]).Destroy;

  for i := 0 to FFruit.Count - 1 do
    for j := 0 to TList(FFruit[i]).Count - 1 do
      TPoint3d(TList(FFruit[i])[j]).Destroy;

  for i := 0 to FSubStems.Count - 1 do
    TStem(FSubStems[i]).Destroy;
  for i := 0 to FStemClones.Count - 1 do
    TStem(FStemClones[i]).Destroy;
  FSections.Destroy;
  FSectionFrameList.Destroy;
  FSubStems.Destroy;
  FStemClones.Destroy;
  FLeafs.Destroy;
  FFruit.Destroy;
  inherited Destroy;
end;

procedure TStem.CreateStructure(const AParent: TStem;
  const ALength, AOffsetChild: real; const Level: Integer);
{ Pre: True
  Post: The complete tree structure has been created. The total number of
  stems and leaves in the tree is known. The length of the stem
  is needed because the amount of substems that will be spawned from
  the current stem depends on the length of the current stem. }
var
  i: Integer; // Loop variables
  SubStemLength: real; // Length of a sub stem
  TotalSubStems: Integer; // The amount of sub stems the current stem
  // will spawn
  SubStem: TStem; // The stem object used to hold spawned sub stems
  FracPos: real; // Holds the current fractional y position along the
  // stem (used when spawning sub stems)
  OffsetSubStem: real; // Holds the current y position along the stem
  // in global coordinates of the sub stem
  // that is about to be spawned
  { StemClones: TStem;
    TotalStemClones: integer;
    OffsetClone,
    CloneLength: real; }
begin
  FOffsetChild := AOffsetChild;
  FParent := AParent;
  FLength := ALength;
  FTree.TotalStems := FTree.TotalStems + 1;

  with FTree.FModelParameters do
  begin
    { Calculate the base length of the current stem
      (is always 0 except for level  0) }
    if Level = 0 then
      FBaseLength := BaseSize * FTree.ScaleTree
    else
      FBaseLength := 0;

    { Spawn sub stems or leaves }

    if (Levelers > Level) then
    begin
      if (Levelers = Level + 1) then
      begin
        { Count the total number of leaves (last level in the recursion) }
        if Level <> 0 then
          FTree.TotalGrownLeaves := FTree.TotalGrownLeaves +
            Round((Round(Leaves * (Health / 100))) * ShapeRatio(4,
            FOffsetChild / FParent.Length) * LeafQuality)
        else
          FTree.TotalGrownLeaves := 0;

        if Level <> 0 then
          FTree.TotalFruit := FTree.TotalFruit +
            Round((Round(Fruits * (Health / 100))) * ShapeRatio(4,
            FOffsetChild / FParent.Length) * LeafQuality)
        else
          FTree.TotalFruit := 0;
      end
      else if (Levelers <> Level + 1) then
      begin { CLONES then SUB STEMS }
        (*
          { Calculate the amount of sub stems the current stem will spawn }
          {Depend on Split or Spawn  Level .. Base Splits }
          { StemClones  TotalStemClones   OffsetClone  CloneLength}
          if Level = 0 then TotalStemClones :=BaseSplits{ Round((1 - BaseSize) * nBranches[1])}
          else  TotalStemClones := Round(nSegSplits[Level] * (1.0 - 0.5 * FOffsetChild / FParent.Length));

          { Add Clones coming out of the current stem }
          for i := 0 to TotalStemClones - 1 do
          begin
          if Level = 0 then FracPos :=0{ BaseSize + ((i + 1) * (1 - BaseSize) / (TotalSubStems + 1))}
          else FracPos := (i + 1) * (1 / (TotalStemClones + 1));

          { Calculate the y position of the sub stem measured from the base
          of the stem (in the global coordinate system) }
          OffsetClone := FracPos * FLength;

          { Calculate the length of the sub stem }
          FLengthCloneMax := nLength[Level ] + GetRandomValue(nLengthV[Level ]);

          if Level = 0 then
          begin
          If (Shape<8) then
          CloneLength := FTree.TrunkLength * FLengthCloneMax *
          ShapeRatio(Shape, (FTree.TrunkLength - OffsetClone)
          / (FTree.TrunkLength - FBaseLength))
          else {Do Prune}
          {r := sqrt ( x2 + y2 + z2 )
          ratio := ( scaletree - z ) / ( scaletree *(1-BaseSize ) )
          inside := [ r / scaletree < PruneWidth *
          ShapeRatio( 8 , ratio ) ]}
          CloneLength := FTree.TrunkLength * FLengthCloneMax *
          ShapeRatio(Shape, (FTree.TrunkLength - OffsetClone)
          / (FTree.TrunkLength - FBaseLength));
          {StrtoFloat(PruneRatioFloatEdit.Text)  0=no use..1.0 :TOTAL USE
          StrtoFloat(ScaleFloatEdit.Text)
          StrtoFloat(ScaleVFloatEdit.Text)
          StrtoFloat(BaseSizeFloatEdit.Text)
          StrtoFloat(Length0FloatEdit.Text)
          StrtoFloat(LengthV0FloatEdit.Text)}
          {
          r := sqrt ( x2 + y2 + z2 )
          ratio := ( scaletree - z ) / ( scaletree *(1-BaseSize ) )
          inside := [ r / scaletree < PruneWidth *ShapeRatio( 8 , ratio ) ]
          where ShapeRatio( 8, ratio ) is defined as
          [ ratio / ( 1 - PruneWidthPeak ) ]PrunePowerHigh
          when ratio < 1 - PruneWidthPeak, or
          [ ( 1 - ratio ) / ( 1 - PruneWidthPeak ) ]PrunePowerLow
          }

          end
          else
          CloneLength := FLengthCloneMax * (FLength - 0.6 * OffsetClone);

          { Spawn the sub stem,
          but only if the sub stem radius is greater than zero }

          StemClones := TStem.Create(FTree);
          FStemClones.Add(StemClones);
          StemClones.CreateStructure(Self, CloneLength, OffsetClone, Level + 1);
          end;
        *)

        { Calculate the amount of sub stems the current stem will spawn }
        { Depend on Split or Spawn  Level .. Base Splits }
        if Level = 0 then
          TotalSubStems := Round((1 - BaseSize) * nBranches[1])
        else if Level = 1 then
          TotalSubStems :=
            Round(nBranches[2] * (0.2 + 0.8 * (FLength / FParent.Length) /
            FParent.LengthChildMax))
        else
          TotalSubStems :=
            Round(nBranches[Level + 1] * (1.0 - 0.5 * FOffsetChild /
            FParent.Length));

        { Add child stems coming out of the current stem }

        for i := 0 to TotalSubStems - 1 do
        begin
          { FTree.TotalTreeLeavesGrown :=
            FTree.TotalTreeLeavesGrown
            + Round(
            (Round(nLeaves[Level] *(Health/100)))
            * ShapeRatio(4, FOffsetChild / FParent.Length)
            * LeafQuality); }
          { Calculate fractional position along the stem }

          if Level = 0 then
            FracPos := BaseSize +
              ((i + 1) * (1 - BaseSize) / (TotalSubStems + 1))
          else
            FracPos := (i + 1) * (1 / (TotalSubStems + 1));

          { Calculate the y position of the sub stem measured from the base
            of the stem (in the global coordinate system) }

          OffsetSubStem := FracPos * FLength;

          { Calculate the length of the sub stem }
          FLengthChildMax := nLength[Level + 1] +
            GetRandomValue(nLengthV[Level + 1]);

          if Level = 0 then
          begin
            If (Shape < 8) then
              SubStemLength := FTree.TrunkLength * FLengthChildMax *
                ShapeRatio(Shape, (FTree.TrunkLength - OffsetSubStem) /
                (FTree.TrunkLength - FBaseLength))
            else { Do Prune }
              { r := sqrt ( x2 + y2 + z2 )
                ratio := ( scaletree - z ) / ( scaletree *(1-BaseSize ) )
                inside := [ r / scaletree < PruneWidth *
                ShapeRatio( 8 , ratio ) ] }
              SubStemLength := FTree.TrunkLength * FLengthChildMax *
                ShapeRatio(Shape, (FTree.TrunkLength - OffsetSubStem) /
                (FTree.TrunkLength - FBaseLength));
            { StrtoFloat(PruneRatioFloatEdit.Text)  0=no use..1.0 :TOTAL USE
              StrtoFloat(ScaleFloatEdit.Text)
              StrtoFloat(ScaleVFloatEdit.Text)
              StrtoFloat(BaseSizeFloatEdit.Text)
              StrtoFloat(Length0FloatEdit.Text)
              StrtoFloat(LengthV0FloatEdit.Text) }
            {
              r := sqrt ( x2 + y2 + z2 )
              ratio := ( scaletree - z ) / ( scaletree *(1-BaseSize ) )
              inside := [ r / scaletree < PruneWidth *ShapeRatio( 8 , ratio ) ]
              where ShapeRatio( 8, ratio ) is defined as
              [ ratio / ( 1 - PruneWidthPeak ) ]PrunePowerHigh
              when ratio < 1 - PruneWidthPeak, or
              [ ( 1 - ratio ) / ( 1 - PruneWidthPeak ) ]PrunePowerLow
            }

          end
          else
            SubStemLength := FLengthChildMax * (FLength - 0.6 * OffsetSubStem);

          { Spawn the sub stem,
            but only if the sub stem radius is greater than zero }

          SubStem := TStem.Create(FTree);
          FSubStems.Add(SubStem);
          SubStem.CreateStructure(Self, SubStemLength, OffsetSubStem,
            Level + 1);

        end;
      end;
    end;
  end;
end;

procedure TStem.Grow(const AFrame: TCoordinateFrame; const ARadius: real;
  const Level: Integer);
{ Pre: AFrame <> nil
  Post: A stem at the recursion level specified by ALevel has been created,
  with length ALength and maximum radius ARadius, using
  coordinate system AFrame. }
var
  i: Integer; // Loop variables
  LocalSectionOrigin: TPoint3d; // Origin of the frame of the section that is
  // being calculated (in local coordinates)
  SectionOrigin: TPoint3d; // Origin of the frame of the section that is
  // being calculated (in global coordinates)
  SectionONB: TONB; // ONB belonging to the frame of the section
  // that is being calculated
  SectionFrame: TCoordinateFrame;
  // Frame of the section that is currently being
  // calculated. The middle of a section is
  // placed on the origin of the frame, and section
  // is in the plane (u, w)
  SectionRadius: real; // Radius of the section that is being calculated
  StemY: real; // Y position along the stem where the current
  // section is located (measured in the coordinate
  // system of the stem, AFrame)
  Angle: real; // The rotation around the v axis for every section
  LocalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem
  GlobalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem in global coordinates
  CurveRotationMatrix: TGLMatrix; // Matrix used to calculate the rotation
  // of the y axis of the local coordinate
  // system of the segment for the curve
  Stepper: Integer; // my added Trunk Texture fix.. that is broken
  a0, a1, a2, swayx, swayy, bx, by, sway_offsetx, sway_offsety: Double;
begin

  { Update the progress bar }
  If TreePlanted then
    FormGlsViewer.ProgressBar.StepIt
  else If DemoTreePlanted then
    ATreeForm.ProgressBar.StepIt
  else
    GlsProgressForm.ProgressBar.StepIt;
  Application.ProcessMessages;

  if GlsUserCancelled then
    Exit;

  FRadius := ARadius;
  StandardOrigin := TPoint3d.Create(0, 0, 0);

  with FTree.FModelParameters do
  begin

    { Now for the amount of sections specified, create a sectionframe and
      create and initialize the sections. }

    SectionFrame := TCoordinateFrame.Create;
    SectionFrame.Assign(AFrame);
    FSectionFrameList.Add(SectionFrame);

    { Create a stem }

    StemY := 0;

    { Calculate the radius of the section. }

    SectionRadius := CalculateSectionRadius(Level, StemY, FLength, FRadius);

    { Create the points that make up the section. }
    CreateSection(SectionFrame, SectionRadius, nVertices[Level]);

    { In case we are creating the trunk, the first segment of the stem is
      (again) divided into nCurveRes[Level] sections. This is done because
      otherwise the exponentional curve at the base of the trunk wouldn't be
      exponentional but linear. }

    if Level = 0 then
    begin

      { We are creating the trunk,
        so divide the first segment in nCurveRes[0] sections. }

      for i := 1 to FLARE_RESOLUTION - 2 do { 1 to 8 }
      begin
        { Calculate the coordinate system of the current section }

        { Calculate the new origin of the next section using the current origin
          and the v vector of the previous section frame ONB }

        LocalSectionOrigin := StandardOrigin.Translate(UnitY,
          FLength / (nCurveRes[Level] * FLARE_RESOLUTION));
        SectionOrigin := SectionFrame.Local2Canonical(LocalSectionOrigin);

        SectionFrame.Origin.Assign(SectionOrigin);
        SectionFrame.SetTransformationLists;

        { Calculate the radius of the section }

        StemY := StemY + (FLength / (nCurveRes[Level] * FLARE_RESOLUTION));

        if FLength = 0 then
          SectionRadius := 0
        else
          SectionRadius := CalculateSectionRadius(Level, StemY / FLength,
            FLength, FRadius);

        { Create the points that make up the section. }

        CreateSection(SectionFrame, SectionRadius, nVertices[Level]);

        { Clean up dynamically created objects }

        SectionOrigin.Destroy;
        LocalSectionOrigin.Destroy;
      end;

      { Recreate the original section frame }

      SectionFrame.Origin.Assign(AFrame.Origin);
      SectionFrame.SetTransformationLists;

    end;

    { Now, for the rest of the sections of the stem, create them and add
      them to the sections list of the stem. }
    { Trying to fix that blip in the Trunk.. dunno how to- this crashes }
    { if Level = 0 then
      begin
      Stepper:=1;
      end else }
    begin
      StemY := 0;
      Stepper := 0;
    end;

    for i := Stepper to nCurveRes[Level] - 1 do
    begin
      { Calculate the coordinate system of the current section }

      { Calculate the new origin of the next section using the current origin
        and the v vector of the previous section frame ONB }

      LocalSectionOrigin := StandardOrigin.Translate(UnitY,
        FLength / nCurveRes[Level]);
      SectionOrigin := SectionFrame.Local2Canonical(LocalSectionOrigin);

      { Calculate the angle over which the y axis of current segment is rotated
        away from the y axis of the previous segment }

      if nCurveBack[Level] <> 0 then
      begin

        { If nCurveBack[Level] is not equal to zero each of the segments in the
          first half of the stem is rotated (nCurve / (nCurveRes / 2)) degrees
          and each in the second half is rotated (nCurveBack / (nCurveRes / 2))
          degrees. }

        if nCurveRes[Level] / (i + 1) < 2 then
          Angle := nCurve[Level] / (nCurveRes[Level] / 2)
        else
          Angle := nCurveBack[Level] / (nCurveRes[Level] / 2);
      end
      else
        Angle := nCurve[Level] / nCurveRes[Level];

      { Blow wind up the Skirts }
      If (WindAnimationActive) then
      begin
        { the sway angles swayx and swayy at unit position Z
          from 0 to 1 of a segment along the length of a stem
          are computed at any "time" (in seconds) using: }
        { a0 = 4 * lengthstem ( 1 - Z ) / radiusZ    (degrees) }
        a0 := 4 * FLength * (1 - i) / SectionRadius;
        a1 := WindStrength / 50 * a0; { (degrees) }
        a2 := (WindGust * WindSign) / 50 * a0 + a1 / 2; { (degrees) }
        { bx = sway_offsetx + radiusstem / lengthstem * time/15     (radians) }
        bx := Random + (FRadius / FLength) * Anime360Counter / 15;
        { by = sway_offsety + radiusstem / lengthstem * time/15     (radians) }
        by := Random + (FRadius / FLength) * Anime360Counter / 15;
        { starts at 0 so only efect is the Random }
        swayx := (a1 * sin(bx) + a2 * sin(0.7 * bx)) / nCurveRes[Level];
        { (degrees) }
        swayy := (a1 * sin(by) + a2 * sin(0.7 * by)) / nCurveRes[Level];
        { (degrees) }
        { The angles sway_offsetx and sway_offsety
          are randomly selected for each stem in the tree... How? }
        { What to do with the X Sway ? }{ shut up compiler till it works }
        If swayy > swayx then;
      end
      else
        swayy := 0;

      Angle := Angle + swayy + GetRandomValue
        (nCurveV[Level] / nCurveRes[Level]);

      { Angle is in degrees, we need the angle in radians, so here we convert the
        angle }

      Angle := DegToRad(Angle);

      { Now calculated the additional angle (in radians) added because of
        vertical attraction. There is no vertical attraction for the trunk 0
        and the main branches 1. }

      if Level > 1 then
        Angle := Angle + CalculateVerticalAttraction(Level, SectionFrame);

      { Calculate the rotationmatrix used to apply a curve to the stem }

      CurveRotationMatrix := TGLMatrix.Create;
      CurveRotationMatrix.SetRotationMatrix(Angle, aX);

      { Calculate the new y vector for the next section }

      LocalY := CurveRotationMatrix.Multiply(UnitY);
      GlobalY := SectionFrame.ONB.LCTransformMatrix.Multiply(LocalY);

      { Create the ONB for the new section }

      SectionONB := TONB.CreateFromUV(SectionFrame.ONB.u, GlobalY);

      { Create a new section frame using the origin and the ONB that have just
        been calculated, and subsequently add the newly created frame to the
        section frame list }

      SectionFrame := TCoordinateFrame.Create(SectionOrigin, SectionONB);
      FSectionFrameList.Add(SectionFrame);

      { Clean up dynamically created temporary objects }

      LocalY.Destroy;
      GlobalY.Destroy;
      LocalSectionOrigin.Destroy;
      SectionOrigin.Destroy;
      SectionONB.Destroy;
      CurveRotationMatrix.Destroy;

      { Calculate the radius of the section }

      StemY := StemY + (FLength / nCurveRes[Level]);

      if FLength = 0 then
        SectionRadius := 0
      else
        SectionRadius := CalculateSectionRadius(Level, StemY / FLength,
          FLength, FRadius);

      { Create the points that make up the section. }

      CreateSection(SectionFrame, SectionRadius, nVertices[Level]);
    end;

    { Spawn sub stems or leaves }

    if Levelers > Level then
    begin

      { We'll need to spawn sub stems or leaves }

      if (Levelers = Level + 1) then
      begin
        if FTree.FModelParameters.nLeaves[Level] > 0 then
          GrowLeaves(Level);
        { Create leaves (last level in the recursion) }
        if FTree.TotalGrownLeaves <> 0 then
          GrowLeaves(Level);
        if FTree.TotalFruit <> 0 then
          GrowFruit(Level);
      end
      else if (Levelers <> Level + 1) then
      begin
        if FTree.FModelParameters.nLeaves[Level] > 0 then
          GrowLeaves(Level);
        { GrowStemCellClones(Level); }
        { Create sub stems }
        GrowSubStems(Level);
      end;
    end;
  end;
  { Clean up dynamically created objects }
  for i := 0 to FSectionFrameList.Count - 1 do
    TCoordinateFrame(FSectionFrameList[i]).Destroy;
end;

{ Calculate the Stem Clone Distribution }
function TStem.CalculateBranchDistribution(const Level: Integer;
  const y, StemLength: real): real;
Begin
  Result := 0.0;
End;

function TStem.CalculateVerticalAttraction(const Level: Integer;
  const SectionFrame: TCoordinateFrame): real;
{ Pre: ALevel > 1 (there is no vertical attraction for the trunk and main branches)
  Ret: An angle in radians that is added to the segments curve angle to simulate
  vertical attraction. }
var
  TransformY, TransformZ: TVector3d; // Vectors holding the unit y and
  // z vectors passed through the
  // current frame transformation
  Declination, Orientation: real; // Declination and orientation
  // of the current segment
begin
  TransformY := SectionFrame.ONB.LCTransformMatrix.Multiply(UnitY);
  TransformZ := SectionFrame.ONB.LCTransformMatrix.Multiply(UnitZ);
  Declination := ArcCos(TransformY.y);
  Orientation := ArcCos(TransformZ.y);
  with FTree.FModelParameters do
    Result := AttractionUp * Declination * Cos(Orientation) / nCurveRes[Level];
  TransformY.Destroy;
  TransformZ.Destroy;
end;

function TStem.ShapeRatio(const Shape: Integer; const Ratio: real): real;
{ Pre: 0 <= Shape <= 8 & 0 <= Ratio <= 1
  Ret: A certain predefined ratio depending on the Shape and Ratio parameter.
  For precise details see "The Creation and Rendering of Realistic Trees". }
begin
  Result := 0;
  case Shape of
    0:
      Result := 0.2 + 0.8 * Ratio;
    1:
      Result := 0.2 + 0.8 * sin(Pi * Ratio);
    2:
      Result := 0.2 + 0.8 * sin(0.5 * Pi * Ratio);
    3:
      Result := 1.0;
    4:
      Result := 0.5 + 0.5 * Ratio;
    5:
      begin
        if Ratio <= 0.7 then
          Result := Ratio / 0.7
        else
          Result := (1.0 - Ratio) / 0.3;
      end;
    6:
      Result := 1.0 - 0.8 * Ratio;
    7:
      begin
        if Ratio <= 0.7 then
          Result := 0.5 + 0.5 * Ratio / 0.7
        else
          Result := 0.5 + 0.5 * (1.0 - Ratio) / 0.3;
      end;
    8 { Result := 1;
      TODO: Use pruning envelope for ShapeRatio(8, Ratio); }
      :
      begin { FTree.FModelParameters. }
        if Ratio < (1 - FTree.FModelParameters.PruneWidthPeak) then
          Result := (Ratio / (1 - FTree.FModelParameters.PruneWidthPeak)) *
            FTree.FModelParameters.PrunePowerHigh
        else
          Result := ((1 - Ratio) / (1 - FTree.FModelParameters.PruneWidthPeak))
            * FTree.FModelParameters.PrunePowerLow;
      end;

    { StrtoFloat(PruneRatioFloatEdit.Text)  0=no use..1.0 :TOTAL USE
      StrtoFloat(ScaleFloatEdit.Text)
      StrtoFloat(ScaleVFloatEdit.Text)
      StrtoFloat(BaseSizeFloatEdit.Text)
      StrtoFloat(Length0FloatEdit.Text)
      StrtoFloat(LengthV0FloatEdit.Text) }
    (*
      r := sqrt ( x2 + y2 + z2 )
      ratio := ( scaletree - z ) / ( scaletree *(1-BaseSize ) )
      inside := [ r / scaletree < PruneWidth *ShapeRatio( 8 {(envelope)}, ratio ) ]
      where ShapeRatio( 8, ratio ) is defined as
      [ ratio / ( 1 - PruneWidthPeak ) ]PrunePowerHigh
      when ratio < 1 - PruneWidthPeak, or
      [ ( 1 - ratio ) / ( 1 - PruneWidthPeak ) ]PrunePowerLow
    *)
  end;
end;

function TStem.CalculateSectionRadius(const Level: Integer;
  const y, StemLength, StemRadius: real): real;
{ Pre: 0 <= y <= 1
  Ret: The radius of the stem at the (normalized) y position along the stem.
  See for the exact calculations the paper "Creation and Rendering of
  Realistic Trees", by Jason Weber and Joseph Penn. }
var
  y2, y3: real; // Temporary variable
  Depth: real; // Scaling factor used for periodic tapering
  TaperY: real; // Tapered radius along at the (normalized) y
  // position along the stem
  UnitTaper: real; // UnitTaper is used to determine the radius of
  // the stem along a specified (normalized)
  // position along the stem.
begin
  with FTree.FModelParameters do
  begin
    { Calculate UnitTaper, a variable used to determine the radius of the
      stem along a specified (normalized) position Z along the stem }

    UnitTaper := 0;

    if (nTaper[Level] >= 0) and (nTaper[Level] < 1) then
      UnitTaper := nTaper[Level]
    else if (nTaper[Level] >= 1) and (nTaper[Level] < 2) then
      UnitTaper := 2 - nTaper[Level]
    else if (nTaper[Level] >= 2) and (nTaper[Level] < 3) then
      UnitTaper := 0;

    TaperY := StemRadius * (1 - (UnitTaper * y));

    if (nTaper[Level] >= 0) and (nTaper[Level] < 1) then
    begin
      Result := TaperY;
    end
    else
    begin
      { (nTaper[ALevel] >= 1) and (nTaper[ALevel] <= 3) }

      { Initialize y2 }

      y2 := (1 - y) * StemLength;

      { Initialize Depth }

      if (nTaper[Level] < 2) or (y2 < TaperY) then
        Depth := 1
      else
        Depth := nTaper[Level] - 2;

      { Initialize y3 }

      if (nTaper[Level] < 2) then
        y3 := y2
      else
        y3 := Abs(y2 - 2 * TaperY * Trunc(y2 / (2 * TaperY) + 0.5));

      { Return the radius }

      if (nTaper[Level] < 2) and (y3 >= TaperY) then
        Result := TaperY
      else
        Result := (1 - Depth) * TaperY + Depth *
          Sqrt(Abs(TaperY * TaperY - (y3 - TaperY) * (y3 - TaperY)));

    end;

    { Calculate flaring }

    if Level = 0 then
    begin
      y2 := 1 - 8 * y;
      if y2 < 0 then
        y2 := 0;
      Result := Result * (Flare * (Power(100, y2) - 1) / 100 + 1);
    end;

  end;
end;

procedure TStem.GrowSubStems(const Level: Integer);
{ Pre: Level >= 0
  Post: The sub stems for the current stem have been created. }
var
  i: Integer; // Loop variables
  LocalSubStemOrigin: TPoint3d; // Origin of the stem that is being spawned
  // (in local coordinates)
  SubStemOrigin: TPoint3d; // Origin of the stem that is being spawned
  // (in global coordinates)
  SubStemONB: TONB; // ONB belonging to the frame of the stem that
  // is being spawned
  SubStemFrame: TCoordinateFrame;
  // Frame of a child stem that is spawned from the
  // current stem. The first section of the child
  // stem is placed on the origin of the frame, and
  // the first section of the child stem in in the
  // plane (u, w), with u parallel to the ground
  // surface, and w pointing to the sky (for a
  // child stem perpendicular to the parent stem).
  FracPos: real; // Holds the current fractional y position along the
  // stem (used when spawning sub stems)
  LocalPos: real; // Holds the current y position along a segment
  CurrentSegment: Integer; // Holds the segment where we reside along the
  // stem (used when spawning sub stems)
  LocalX, LocalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem
  GlobalX, GlobalY: TVector3d;
  // Vectors holding the x and y vector respectively
  // of a local coordinatesystem in global coordinates
  SubStemLength: real; // Length of a sub stem
  SubStemRadius: real; // Radius of a sub stem
  TotalSubStems: Integer; // The amount of sub stems the current stem
  // will spawn
  StemDownAngle: real; // The angle between the current stem and the
  // sub stem
  StemRotateAngle: real; // The angle between the previous sub stem
  // and the currently being spawned sub stem
  LocalRadius: real; // The radius of the stem at the spawning position
  // of the substem
  DownAngleRotationMatrix: TGLMatrix; // Matrix used to calculate the down angle
  // rotation of the y axis of a sub stem
  RotateRotationMatrix: TGLMatrix; // Matrix used to calculate the nRotate
  // rotation of the x axis of a sub stem
  // around the y axis of the parent stem
  OffsetSubStem: real; // Holds the current y position along the stem
  // in global coordinates of the sub stem
  // that is about to be spawned
  TempONB: TONB;
begin
  with FTree.ModelParameters do
  begin
    { Initialize the angle about the y axis of the parent relative to the
      previous sub stem }

    StemRotateAngle := Random * 2 * Pi;

    { Calculate the amount of sub stems the current stem will spawn }

    TotalSubStems := FSubStems.Count;

    { Add child stems coming out of the current stem }

    for i := 0 to TotalSubStems - 1 do
    begin

      { Calculate fractional position along the stem }

      if Level = 0 then
        FracPos := BaseSize + ((i + 1) * (1 - BaseSize) / (TotalSubStems + 1))
      else
        FracPos := (i + 1) * (1 / (TotalSubStems + 1));

      { Calculate the y position of the sub stem measured from the base
        of the stem (in the global coordinate system) }

      OffsetSubStem := FracPos * FLength;

      { Calculate in which segment we reside }

      CurrentSegment := Trunc(FracPos / (1 / nCurveRes[Level])) + 1;

      { Calculate the position within the segment (in meters).
        CurrentSegment * (1 / nCurveRes[Level]) delivers us the fractional y
        position of the origin of the current segment along the stem.
        Multiplying this value with FLength results in the global y position
        in meters of the segment. Subtracting this value from the current
        global y position results in the position within the segment. }

      LocalPos := OffsetSubStem - (CurrentSegment - 1) *
        (FLength / nCurveRes[Level]);

      { Create the origin point of the sub stem frame. The origin is calculated
        from the origin of the local frame of the current segment and the
        local y position within the current segment. }

      LocalSubStemOrigin := StandardOrigin.Translate(UnitY, LocalPos);
      if CurrentSegment = 1 then
        SubStemOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalSubStemOrigin)
      else
        SubStemOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalSubStemOrigin);

      { Create the ONB of the sub stem frame. The x vector is pointing towards
        the watcher, and the y vector is pointing upwards. }

      { Create a rotationmatrix that will be used to rotate the sub stem
        frame around the x-axis }

      DownAngleRotationMatrix := TGLMatrix.Create;

      { Calculate the angle between the current stem and the sub stem and
        use this angle to create a rotation matrix }
      if nDownAngleV[Level + 1] >= 0 then
        StemDownAngle := nDownAngle[Level + 1] +
          GetRandomValue(nDownAngleV[Level + 1])
      else { nDownAngleV[1] < 0 }
        StemDownAngle := nDownAngle[Level + 1] + nDownAngleV[Level + 1] *
          (1 - 2 * ShapeRatio(0, (FLength - OffsetSubStem) /
          (FLength - FBaseLength)));

      { Convert to radians }

      StemDownAngle := DegToRad(StemDownAngle);

      DownAngleRotationMatrix.SetRotationMatrix(StemDownAngle, aX);

      { Now calculate the rotation angle around the y axis of the parent
        relative to the previous sub stem }

      RotateRotationMatrix := TGLMatrix.Create;

      if nRotate[Level + 1] >= 0 then
        StemRotateAngle := StemRotateAngle +
          DegToRad(nRotate[Level + 1] + GetRandomValue(nRotateV[Level + 1]))
      else
        StemRotateAngle := StemRotateAngle +
          DegToRad(180 + nRotate[Level + 1] +
          GetRandomValue(nRotateV[Level + 1]));

      RotateRotationMatrix.SetRotationMatrix(StemRotateAngle, aY);

      { Calculate the new x and y vector for the sub stem, using the
        transformation list that has just been calculated }

      { LocalX now contains the (local) x vector of the frame we are about to
        create. We still have to rotate this x vector around the y axis of the
        current (parent) stem, using the previously created rotation matrix }

      LocalX := RotateRotationMatrix.Multiply(UnitX);

      { Temp now contains the rotated (local) x vector. We'll now convert
        this vector to global coordinates. }

      GlobalX := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
        .ONB.LCTransformMatrix.Multiply(LocalX));

      TempONB := TONB.CreateFromUV(GlobalX,
        TCoordinateFrame(FSectionFrameList[CurrentSegment]).ONB.v);

      LocalY := DownAngleRotationMatrix.Multiply(UnitY);

      { LocalY now contains the y axis of the sub stem local coordinate system
        adjusted for the rotate angle (in local coordinates) }

      GlobalY := (TempONB.LCTransformMatrix.Multiply(LocalY));

      SubStemONB := TONB.CreateFromUV(GlobalX, GlobalY);

      { Create the sub stem frame }

      SubStemFrame := TCoordinateFrame.Create(SubStemOrigin, SubStemONB);

      { Calculate the length of the sub stem }
      FLengthChildMax := nLength[Level + 1] +
        GetRandomValue(nLengthV[Level + 1]);

      if Level = 0 then
        SubStemLength := FTree.TrunkLength * FLengthChildMax *
          ShapeRatio(Shape, (FTree.TrunkLength - OffsetSubStem) /
          (FTree.TrunkLength - FBaseLength))
      else
        SubStemLength := FLengthChildMax * (FLength - 0.6 * OffsetSubStem);

      { Calculate the radius of the sub stem }

      if FLength = 0 then
        SubStemRadius := 0
      else
        SubStemRadius := FRadius * Power((SubStemLength / FLength), RatioPower);

      { Check if the calculated radius is greater than the radius of the current
        stem at the spawning position }

      LocalRadius := CalculateSectionRadius(Level, FracPos, FLength, FRadius);

      if (SubStemRadius > LocalRadius) or (SubStemRadius = 0) then
        SubStemRadius := LocalRadius;

      { Spawn the sub stem, but only if the sub stem radius is greater than
        zero }

      TStem(FSubStems[i]).Grow(SubStemFrame, SubStemRadius, Level + 1);

      { Clean up dynamically created temporary variables }

      LocalX.Destroy;
      LocalY.Destroy;
      GlobalX.Destroy;
      GlobalY.Destroy;
      TempONB.Destroy;
      LocalSubStemOrigin.Destroy;
      SubStemOrigin.Destroy;
      SubStemONB.Destroy;
      DownAngleRotationMatrix.Destroy;
      RotateRotationMatrix.Destroy;
      SubStemFrame.Destroy;

    end;
  end;
end;

procedure TStem.GrowStemCellClones(const Level: Integer);
{ Pre: Level >= 0
  Post: The Clones for the current stem have been created. }
var
  i: Integer; // Loop variables
  LocalCloneOrigin: TPoint3d; // Origin of the stem that is being spawned
  // (in local coordinates)
  CloneOrigin: TPoint3d; // Origin of the stem that is being spawned
  // (in global coordinates)
  CloneONB: TONB; // ONB belonging to the frame of the stem that
  // is being spawned
  CloneFrame: TCoordinateFrame;
  // Frame of a child stem that is spawned from the
  // current stem. The first section of the child
  // stem is placed on the origin of the frame, and
  // the first section of the child stem in in the
  // plane (u, w), with u parallel to the ground
  // surface, and w pointing to the sky (for a
  // child stem perpendicular to the parent stem).
  FracPos: real; // Holds the current fractional y position along the
  // stem (used when spawning sub stems)
  LocalPos: real; // Holds the current y position along a segment
  CurrentSegment: Integer; // Holds the segment where we reside along the
  // stem (used when spawning sub stems)
  LocalX, LocalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem
  GlobalX, GlobalY: TVector3d;
  // Vectors holding the x and y vector respectively
  // of a local coordinatesystem in global coordinates
  CloneLength: real; // Length of a sub stem
  CloneRadius: real; // Radius of a sub stem
  TotalClones: Integer; // The amount of Clones the current stem
  // will spawn
  StemDownAngle: real; // The angle between the current stem and the
  // sub stem
  StemRotateAngle: real; // The angle between the previous sub stem
  // and the currently being spawned sub stem
  LocalRadius: real; // The radius of the stem at the spawning position
  // of the substem
  DownAngleRotationMatrix: TGLMatrix; // Matrix used to calculate the down angle
  // rotation of the y axis of a sub stem
  RotateRotationMatrix: TGLMatrix; // Matrix used to calculate the nRotate
  // rotation of the x axis of a sub stem
  // around the y axis of the parent stem
  OffsetClone: real; // Holds the current y position along the stem
  // in global coordinates of the sub stem
  // that is about to be spawned
  TempONB: TONB;
begin
  with FTree.ModelParameters do
  begin
    { Initialize the angle about the y axis of the parent relative to the
      previous sub stem }
    { If Level=0 then
      begin
      For i:=0 to BaseSplits do begin
      if BaseSplits > 1 then        BaseSplitsCounter
      StemRotateAngle :=( (i*360)/(BaseSplits+1) else
      end; }

    StemRotateAngle := Random * 2 * Pi;

    { Calculate the amount of Clones the current stem will spawn }

    TotalClones := FStemClones.Count;

    { Add child stems coming out of the current stem }

    for i := 0 to TotalClones - 1 do
    begin

      { Calculate fractional position along the stem }

      if Level = 0 then
        FracPos :=
          0 { BaseSize + ((i + 1) * (1 - BaseSize) / (TotalClones + 1)) }
      else
        FracPos := (i + 1) * (1 / (TotalClones + 1));

      { Calculate the y position of the sub stem measured from the base
        of the stem (in the global coordinate system) }

      OffsetClone := FracPos * FLength;

      { Calculate in which segment we reside }

      CurrentSegment := Trunc(FracPos / (1 / nCurveRes[Level])) + 1;

      { Calculate the position within the segment (in meters).
        CurrentSegment * (1 / nCurveRes[Level]) delivers us the fractional y
        position of the origin of the current segment along the stem.
        Multiplying this value with FLength results in the global y position
        in meters of the segment. Subtracting this value from the current
        global y position results in the position within the segment. }

      LocalPos := OffsetClone - (CurrentSegment - 1) *
        (FLength / nCurveRes[Level]);

      { Create the origin point of the sub stem frame. The origin is calculated
        from the origin of the local frame of the current segment and the
        local y position within the current segment. }

      LocalCloneOrigin := StandardOrigin.Translate(UnitY, LocalPos);
      if CurrentSegment = 1 then
        CloneOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalCloneOrigin)
      else
        CloneOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalCloneOrigin);

      { Create the ONB of the sub stem frame. The x vector is pointing towards
        the watcher, and the y vector is pointing upwards. }

      { Create a rotationmatrix that will be used to rotate the sub stem
        frame around the x-axis }

      DownAngleRotationMatrix := TGLMatrix.Create;

      { Calculate the angle between the current stem and the sub stem and
        use this angle to create a rotation matrix }
      if nDownAngleV[Level] >= 0 then
        StemDownAngle := nDownAngle[Level] + GetRandomValue(nDownAngleV[Level])
      else { nDownAngleV[1] < 0 }
        StemDownAngle := nDownAngle[Level] + nDownAngleV[Level] *
          (1 - 2 * ShapeRatio(0, (FLength - OffsetClone) /
          (FLength - FBaseLength)));

      { Convert to radians }

      StemDownAngle := DegToRad(StemDownAngle);

      DownAngleRotationMatrix.SetRotationMatrix(StemDownAngle, aX);

      { Now calculate the rotation angle around the y axis of the parent
        relative to the previous sub stem }

      RotateRotationMatrix := TGLMatrix.Create;

      if nRotate[Level] >= 0 then
        StemRotateAngle := StemRotateAngle +
          DegToRad(nRotate[Level] + GetRandomValue(nRotateV[Level]))
      else
        StemRotateAngle := StemRotateAngle +
          DegToRad(180 + nRotate[Level] + GetRandomValue(nRotateV[Level]));

      RotateRotationMatrix.SetRotationMatrix(StemRotateAngle, aY);

      { Calculate the new x and y vector for the sub stem, using the
        transformation list that has just been calculated }

      { LocalX now contains the (local) x vector of the frame we are about to
        create. We still have to rotate this x vector around the y axis of the
        current (parent) stem, using the previously created rotation matrix }

      LocalX := RotateRotationMatrix.Multiply(UnitX);

      { Temp now contains the rotated (local) x vector. We'll now convert
        this vector to global coordinates. }

      GlobalX := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
        .ONB.LCTransformMatrix.Multiply(LocalX));

      TempONB := TONB.CreateFromUV(GlobalX,
        TCoordinateFrame(FSectionFrameList[CurrentSegment]).ONB.v);

      LocalY := DownAngleRotationMatrix.Multiply(UnitY);

      { LocalY now contains the y axis of the sub stem local coordinate system
        adjusted for the rotate angle (in local coordinates) }

      GlobalY := (TempONB.LCTransformMatrix.Multiply(LocalY));

      CloneONB := TONB.CreateFromUV(GlobalX, GlobalY);

      { Create the sub stem frame }

      CloneFrame := TCoordinateFrame.Create(CloneOrigin, CloneONB);

      { Calculate the length of the sub stem }
      FLengthCloneMax := nLength[Level] + GetRandomValue(nLengthV[Level]);

      if Level = 0 then
        CloneLength := FTree.TrunkLength * FLengthCloneMax *
          ShapeRatio(Shape, (FTree.TrunkLength - OffsetClone) /
          (FTree.TrunkLength - FBaseLength))
      else
        CloneLength := FLengthCloneMax * (FLength - 0.6 * OffsetClone);

      { Calculate the radius of the sub stem }

      if FLength = 0 then
        CloneRadius := 0
      else
        CloneRadius := FRadius * Power((CloneLength / FLength), RatioPower);

      { Check if the calculated radius is greater than the radius of the current
        stem at the spawning position }

      LocalRadius := CalculateSectionRadius(Level, FracPos, FLength, FRadius);

      if (CloneRadius > LocalRadius) or (CloneRadius = 0) then
        CloneRadius := LocalRadius;

      { Spawn the sub stem, but only if the sub stem radius is greater than
        zero }

      TStem(StemClones[i]).Grow(CloneFrame, CloneRadius, Level);

      { Clean up dynamically created temporary variables }

      LocalX.Destroy;
      LocalY.Destroy;
      GlobalX.Destroy;
      GlobalY.Destroy;
      TempONB.Destroy;
      LocalCloneOrigin.Destroy;
      CloneOrigin.Destroy;
      CloneONB.Destroy;
      DownAngleRotationMatrix.Destroy;
      RotateRotationMatrix.Destroy;
      CloneFrame.Destroy;
    end;
  end;
end;

procedure TStem.GrowLeaves(const Level: Integer);
{ Pre: ALevel >= 0
  Post: The leaves for the current stem have been created. }
var
  LocalLeafOrigin: TPoint3d; // Origin of the leaf that is being spawned
  // (in local coordinates)
  LeafOrigin: TPoint3d; // Origin of the leaf that is being spawned
  // (in global coordinates)
  LeafONB: TONB; // ONB belonging to the frame of the leaf that
  // is being spawned
  LeafFrame: TCoordinateFrame; // Frame of a leaf that is spawned from the
  // current stem.
  TotalLocalLeaves: Integer; // The amount of leaves the current stem will
  // spawn
  i: Integer; // Loop variables
  FracPos: real; // Holds the current fractional y position along the
  // stem (used when spawning sub stems)
  CurrentSegment: Integer; // Holds the segment where we reside along the
  // stem (used when spawning sub stems)
  LocalPos: real; // Holds the current y position along a segment
  RotateRotationMatrix: TGLMatrix; // Matrix used to calculate the nRotate
  // rotation of the x axis of a sub stem
  // around the y axis of the parent stem
  DownAngleRotationMatrix: TGLMatrix; // Matrix used to calculate the down angle
  // rotation of the y axis of a sub stem
  LeafRotateAngle: real; // The angle between the previous leaf
  // and the currently being spawned leaf
  LeafDownAngle: real; // The angle between the current stem and the
  // leaf
  LocalX, LocalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem
  GlobalX, GlobalY: TVector3d;
  // Vectors holding the x and y vector respectively
  // of a local coordinatesystem in global coordinates
  Temp: TVector3d; // Vector used to store temporary results of
  // bigger calculations
  OffsetLeaf: real;
begin
  with FTree.ModelParameters do
  begin
    TotalLocalLeaves := 0;
    if FTree.FModelParameters.nLeaves[Level] > 0 then
      TotalLocalLeaves := TotalLocalLeaves + nLeaves[Level]
      { + Round(
        (Round(nLeaves[Level] *(Health/100)))
        * ShapeRatio(4, FOffsetChild / FParent.Length)
        * LeafQuality) };
    if (Levelers = Level + 1) then
      TotalLocalLeaves := TotalLocalLeaves +
        Round((Round(Leaves * (Health / 100))) * ShapeRatio(4,
        FOffsetChild / FParent.Length) * LeafQuality);

    LeafRotateAngle := Random * 2 * Pi;
    { TreeLeavesAddedTotal:=TreeLeavesAddedTotal+ TotalLocalLeaves; }
    for i := 0 to TotalLocalLeaves - 1 do
    begin

      { Calculate fractional position along the stem }

      if Level = 0 then
        FracPos := BaseSize +
          ((i + 1) * (1 - BaseSize) / (TotalLocalLeaves + 1))
      else
        FracPos := (i + 1) * (1 / (TotalLocalLeaves + 1));

      { Calculate the y position of the leaf measured from the base
        of the stem (in the global coordinate system) }

      OffsetLeaf := FracPos * FLength;

      { Calculate in which segment we reside }

      CurrentSegment := Trunc(FracPos / (1 / nCurveRes[Level])) + 1;

      { Calculate the position within the segment (in meters).
        CurrentSegment * (1 / nCurveRes[0]) delivers us the fractional y
        position of the origin of the current segment along the stem.
        Multiplying this value with FLength results in the global y position
        in meters of the segment. Subtracting this value from the current
        global y position results in the position within the segment. }

      LocalPos := FracPos * FLength - (CurrentSegment - 1) *
        (FLength / nCurveRes[Level]);

      { Create the leaf frame }

      { Create the origin point of the sub stem frame. The origin is calculated
        from the origin of the local frame of the current segment and the
        local y position within the current segment. }

      LocalLeafOrigin := StandardOrigin.Translate(UnitY, LocalPos);
      if CurrentSegment = 1 then
        LeafOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalLeafOrigin)
      else
        LeafOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalLeafOrigin);

      { Create the ONB of the sub stem frame. The x vector is pointing towards
        the watcher, and the y vector is pointing upwards. }

      { *** Not in the model: when the season is autumn, some leaves will
        be lying on the ground. }

      { Create a rotationmatrix that will be used to rotate the sub stem
        around the y axis of the parent relative to the previous sub stem }

      RotateRotationMatrix := TGLMatrix.Create;

      { Now calculate the rotation angle around the y axis of the parent
        relative to the previous sub stem }

      if nRotate[Level] >= 0 then
        LeafRotateAngle := LeafRotateAngle +
          DegToRad(nRotate[Level] + GetRandomValue(nRotateV[Level]))
      else
        LeafRotateAngle := LeafRotateAngle +
          DegToRad(180 + nRotate[Level] + GetRandomValue(nRotateV[Level]));

      RotateRotationMatrix.SetRotationMatrix(LeafRotateAngle, aY);

      if (FTree.Season = 3) and (Random(3) = 0) then
      begin

        { Adjust the y coordinates of the axes of the local coordinate system
          and of the origin of the local coordinate system }

        LeafFrame := TCoordinateFrame.Create;
        LeafFrame.Origin.Assign(LeafOrigin);

        { Scatter the leaves around over the ground }

        LeafFrame.Origin.y := 0;
        LeafFrame.Origin.x := LeafFrame.Origin.x * (1 + Random * 3);
        LeafFrame.Origin.z := LeafFrame.Origin.z * (1 + Random * 3);

        GlobalX := RotateRotationMatrix.Multiply(UnitX);

        LeafFrame.ONB := TONB.CreateFromVW(GlobalX, UnitY);

        LeafFrame.SetTransformationLists;

        { Clean up dynamically created temporary variables }

        GlobalX.Destroy;
        LocalLeafOrigin.Destroy;
        LeafOrigin.Destroy;
        RotateRotationMatrix.Destroy;
      end
      else
      begin

        { Create a rotationmatrix that will be used to rotate the sub stem
          frame around the x-axis }

        DownAngleRotationMatrix := TGLMatrix.Create;

        { Calculate the angle between the current stem and the sub stem and
          use this angle to create a rotation matrix }

        if nDownAngleV[Level] >= 0 then
          LeafDownAngle := nDownAngle[Level] +
            GetRandomValue(nDownAngleV[Level])
        else { nDownAngleV[1] < 0 }
          LeafDownAngle := nDownAngle[Level] + nDownAngleV[Level] *
            (1 - 2 * ShapeRatio(0, (FLength - OffsetLeaf) /
            (FLength - FBaseLength)));

        { Convert to radians }

        LeafDownAngle := DegToRad(LeafDownAngle);

        DownAngleRotationMatrix.SetRotationMatrix(LeafDownAngle, aZ);

        { Calculate the new x and y vector for the sub stem, using the
          transformation list that has just been calculated }
        { 1 moves them to odd places }
        LocalX := UnitZ.Multiply(-1);

        { LocalX now contains the (local) x vector of the frame we are about to
          create. We still have to rotate this x vector around the y axis of the
          current (parent) stem, using the previously created rotation matrix }

        Temp := RotateRotationMatrix.Multiply(LocalX);

        { Temp now contains the rotated (local) x vector. We'll now convert
          this vector to global coordinates. }

        GlobalX := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
          .ONB.LCTransformMatrix.Multiply(Temp));

        Temp.Destroy;
        Temp := DownAngleRotationMatrix.Multiply(UnitY);

        { Temp now contains the y axis of the sub stem local coordinate system
          adjusted for the down angle (in local coordinates) }

        LocalY := RotateRotationMatrix.Multiply(Temp);

        { LocalY now contains the y axis of the sub stem local coordinate system
          adjusted for the rotate angle (in local coordinates) }

        GlobalY := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
          .ONB.LCTransformMatrix.Multiply(LocalY));

        LeafONB := TONB.CreateFromUV(GlobalX, GlobalY);

        LeafFrame := TCoordinateFrame.Create(LeafOrigin, LeafONB);

        { Clean up dynamically created temporary variables }

        LocalX.Destroy;
        LocalY.Destroy;
        GlobalX.Destroy;
        GlobalY.Destroy;
        Temp.Destroy;
        DownAngleRotationMatrix.Destroy;
        RotateRotationMatrix.Destroy;
        LocalLeafOrigin.Destroy;
        LeafOrigin.Destroy;
        LeafONB.Destroy;
      end;

      { We will now adjust the leaf orientation, so that it is facing upwards
        and outwards, to optimize the available direct sunlight and scattered
        sky light }
      { LeafFrame:=AdjustLeafFrame(LeafFrame); }

      { And finally, create the leaf itself }
      CreateLeaf(LeafShape, LeafFrame);

      { Update the progress bar }
      If TreePlanted then
        FormGlsViewer.ProgressBar.StepIt
      else If DemoTreePlanted then
        ATreeForm.ProgressBar.StepIt
      else
        GlsProgressForm.ProgressBar.StepIt;
      Application.ProcessMessages;
      { Clean up dynamically created temporary variables }
      LeafFrame.Destroy;
    end;
  end;
end;

procedure TStem.GrowFruit(const Level: Integer);
{ Pre: ALevel >= 0
  Post: The leaves for the current stem have been created. }
var
  LocalFruitOrigin: TPoint3d; // Origin of the leaf that is being spawned
  // (in local coordinates)
  FruitOrigin: TPoint3d; // Origin of the leaf that is being spawned
  // (in global coordinates)
  FruitONB: TONB; // ONB belonging to the frame of the leaf that
  // is being spawned
  FruitFrame: TCoordinateFrame; // Frame of a leaf that is spawned from the
  // current stem.
  TotalFruit: Integer; // The amount of leaves the current stem will
  // spawn
  i: Integer; // Loop variables
  FracPos: real; // Holds the current fractional y position along the
  // stem (used when spawning sub stems)
  CurrentSegment: Integer; // Holds the segment where we reside along the
  // stem (used when spawning sub stems)
  LocalPos: real; // Holds the current y position along a segment
  RotateRotationMatrix: TGLMatrix; // Matrix used to calculate the nRotate
  // rotation of the x axis of a sub stem
  // around the y axis of the parent stem
  DownAngleRotationMatrix: TGLMatrix; // Matrix used to calculate the down angle
  // rotation of the y axis of a sub stem
  FruitRotateAngle: real; // The angle between the previous leaf
  // and the currently being spawned leaf
  FruitDownAngle: real; // The angle between the current stem and the
  // leaf
  LocalX, LocalY: TVector3d; // Vectors holding the x and y vector respectively
  // of a local coordinatesystem
  GlobalX, GlobalY: TVector3d;
  // Vectors holding the x and y vector respectively
  // of a local coordinatesystem in global coordinates
  Temp: TVector3d; // Vector used to store temporary results of
  // bigger calculations
  OffsetFruit: real;
begin
  with FTree.ModelParameters do
  begin
    TotalFruit :=
    { Round(FruitPercent * ShapeRatio(4, FOffsetChild / FParent.Length) * LeafQuality); }
      Round((Round(Fruits * (Health / 100))) * ShapeRatio(4,
      FOffsetChild / FParent.Length) * LeafQuality);
    FruitRotateAngle := Random * 2 * Pi;

    for i := 0 to TotalFruit - 1 do
    begin

      { Calculate fractional position along the stem }

      if Level = 0 then
        FracPos := BaseSize + ((i + 1) * (1 - BaseSize) / (TotalFruit + 1))
      else
        FracPos := (i + 1) * (1 / (TotalFruit + 1));

      { Calculate the y position of the leaf measured from the base
        of the stem (in the global coordinate system) }

      OffsetFruit := FracPos * FLength;

      { Calculate in which segment we reside }

      CurrentSegment := Trunc(FracPos / (1 / nCurveRes[Level])) + 1;

      { Calculate the position within the segment (in meters).
        CurrentSegment * (1 / nCurveRes[0]) delivers us the fractional y
        position of the origin of the current segment along the stem.
        Multiplying this value with FLength results in the global y position
        in meters of the segment. Subtracting this value from the current
        global y position results in the position within the segment. }

      LocalPos := FracPos * FLength - (CurrentSegment - 1) *
        (FLength / nCurveRes[Level]);

      { Create the leaf frame }

      { Create the origin point of the sub stem frame. The origin is calculated
        from the origin of the local frame of the current segment and the
        local y position within the current segment. }

      LocalFruitOrigin := StandardOrigin.Translate(UnitY, LocalPos);
      if CurrentSegment = 1 then
        FruitOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalFruitOrigin)
      else
        FruitOrigin := TCoordinateFrame(FSectionFrameList[CurrentSegment - 1])
          .Local2Canonical(LocalFruitOrigin);

      { Create the ONB of the sub stem frame. The x vector is pointing towards
        the watcher, and the y vector is pointing upwards. }

      { *** Not in the model: when the season is autumn, some leaves will
        be lying on the ground. }

      { Create a rotationmatrix that will be used to rotate the sub stem
        around the y axis of the parent relative to the previous sub stem }

      RotateRotationMatrix := TGLMatrix.Create;

      { Now calculate the rotation angle around the y axis of the parent
        relative to the previous sub stem }

      if nRotate[Level] >= 0 then
        FruitRotateAngle := FruitRotateAngle +
          DegToRad(nRotate[Level] + GetRandomValue(nRotateV[Level]))
      else
        FruitRotateAngle := FruitRotateAngle +
          DegToRad(180 + nRotate[Level] + GetRandomValue(nRotateV[Level]));

      RotateRotationMatrix.SetRotationMatrix(FruitRotateAngle, aY);

      if (FTree.Season = 3) and (Random(3) = 0) then
      begin { Make Fruit Fall }

        { Adjust the y coordinates of the axes of the local coordinate system
          and of the origin of the local coordinate system }

        FruitFrame := TCoordinateFrame.Create;
        FruitFrame.Origin.Assign(FruitOrigin);

        { Scatter the Fruit around over the ground }

        FruitFrame.Origin.y := 0;
        FruitFrame.Origin.x := FruitFrame.Origin.x * (1 + Random * 3);
        FruitFrame.Origin.z := FruitFrame.Origin.z * (1 + Random * 3);

        GlobalX := RotateRotationMatrix.Multiply(UnitX);

        FruitFrame.ONB := TONB.CreateFromVW(GlobalX, UnitY);

        FruitFrame.SetTransformationLists;

        { Clean up dynamically created temporary variables }

        GlobalX.Destroy;
        LocalFruitOrigin.Destroy;
        FruitOrigin.Destroy;
        RotateRotationMatrix.Destroy;
      end
      else
      begin

        { Create a rotationmatrix that will be used to rotate the sub stem
          frame around the x-axis }

        DownAngleRotationMatrix := TGLMatrix.Create;

        { Calculate the angle between the current stem and the sub stem and
          use this angle to create a rotation matrix }

        if nDownAngleV[Level] >= 0 then
          FruitDownAngle := nDownAngle[Level] +
            GetRandomValue(nDownAngleV[Level])
        else { nDownAngleV[1] < 0 }
          FruitDownAngle := nDownAngle[Level] + nDownAngleV[Level] *
            (1 - 2 * ShapeRatio(0, (FLength - OffsetFruit) /
            (FLength - FBaseLength)));

        { Convert to radians }

        FruitDownAngle := DegToRad(FruitDownAngle);

        DownAngleRotationMatrix.SetRotationMatrix(FruitDownAngle, aZ);

        { Calculate the new x and y vector for the sub stem, using the
          transformation list that has just been calculated }

        LocalX := UnitZ.Multiply(-1);

        { LocalX now contains the (local) x vector of the frame we are about to
          create. We still have to rotate this x vector around the y axis of the
          current (parent) stem, using the previously created rotation matrix }

        Temp := RotateRotationMatrix.Multiply(LocalX);

        { Temp now contains the rotated (local) x vector. We'll now convert
          this vector to global coordinates. }

        GlobalX := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
          .ONB.LCTransformMatrix.Multiply(Temp));

        Temp.Destroy;
        Temp := DownAngleRotationMatrix.Multiply(UnitY);

        { Temp now contains the y axis of the sub stem local coordinate system
          adjusted for the down angle (in local coordinates) }

        LocalY := RotateRotationMatrix.Multiply(Temp);

        { LocalY now contains the y axis of the sub stem local coordinate system
          adjusted for the rotate angle (in local coordinates) }

        GlobalY := (TCoordinateFrame(FSectionFrameList[CurrentSegment])
          .ONB.LCTransformMatrix.Multiply(LocalY));

        FruitONB := TONB.CreateFromUV(GlobalX, GlobalY);

        FruitFrame := TCoordinateFrame.Create(FruitOrigin, FruitONB);

        { Clean up dynamically created temporary variables }

        LocalX.Destroy;
        LocalY.Destroy;
        GlobalX.Destroy;
        GlobalY.Destroy;
        Temp.Destroy;
        DownAngleRotationMatrix.Destroy;
        RotateRotationMatrix.Destroy;
        LocalFruitOrigin.Destroy;
        FruitOrigin.Destroy;
        FruitONB.Destroy;
      end;

      { We will now adjust the leaf orientation, so that it is facing upwards
        and outwards, to optimize the available direct sunlight and scattered
        sky light }

      { And finally, create the Fruit itself }

      CreateFruit(FruitFrame);

      { Update the progress bar }
      If TreePlanted then
        FormGlsViewer.ProgressBar.StepIt
      else If DemoTreePlanted then
        ATreeForm.ProgressBar.StepIt
      else
        GlsProgressForm.ProgressBar.StepIt;
      Application.ProcessMessages;

      { Clean up dynamically created temporary variables }
      FruitFrame.Destroy;
    end;
  end;
end;

procedure TStem.CreateSection(const AFrame: TCoordinateFrame;
  SectionRadius: real; const Vertices: Integer);
{ Pre: AFrame <> nil
  Post: The FPoints list has been filled with points that make up the section.
  The amount of points depends on the quality set by the user. }
var
  i: Integer; // Loop variable
  Angle: real; // Angle holds the angle in radians
  // between two points in the section.
  LocalPoint: TPoint3d; // LocalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // local coordinates)
  GlobalPoint: TPoint3d; // GlobalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // global coordinates)
  Points: TList; // List of points
  LobedSectionRadius: real; // The sectionradius rescaled using
  // the lobing parameters
begin
  Points := TList.Create;

  { Lame but effective hack to prevent empty triangles }
  if SectionRadius = 0 then
    SectionRadius := 0.0001;

  Angle := (2 * Pi) / Vertices;
  for i := 0 to Vertices - 1 do
  begin

    { Apply lobing to the section }

    with FTree.ModelParameters do
      LobedSectionRadius := SectionRadius *
        (1.0 + LobeDepth * sin(Lobes * i * Angle));

    LocalPoint := TPoint3d.Create;
    LocalPoint.x := Cos(i * Angle) * LobedSectionRadius;
    LocalPoint.y := 0;
    LocalPoint.z := sin(i * Angle) * LobedSectionRadius;

    GlobalPoint := AFrame.Local2Canonical(LocalPoint);
    Points.Add(GlobalPoint);

    { Update the maximum x, y and z values of the tree }

    with FTree do
    begin
      if GlobalPoint.x > MaxX then
        MaxX := GlobalPoint.x;
      if GlobalPoint.y > MaxY then
        MaxY := GlobalPoint.y;
      if GlobalPoint.z > MaxZ then
        MaxZ := GlobalPoint.z;
    end;

    LocalPoint.Destroy;
  end;
  Sections.Add(Points);
end;

procedure TStem.CreateLeaf(const LeafShape: Integer;
  const AFrame: TCoordinateFrame);
var
  i: Integer; // Loop variable
  Points: TList; // List of points
  LocalPoint: TPoint3d; // LocalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // local coordinates)
  GlobalPoint: TPoint3d; // GlobalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // global coordinates)
begin
  Points := TList.Create;
  { Depending on range use shape defined
    or simplest..     to less shapes+bigger leaves }
  { THIS is WHERE the Leaf Patterns would go
    Thus be able to make Shaped Globs of Leaves...
    But as Multiple .. Counts WRONG }
  for i := 0 to LeafPolygonVerticesNumber[LeafShape] do
  begin
    LocalPoint := TPoint3d.Create;
    LocalPoint.x := LeafPolygonVertices[LeafShape, i, 0];
    { Get out of Branch Axis... invisible Stalk }
    LocalPoint.y := LeafPolygonVertices[LeafShape, i, 1] +
      FTree.FModelParameters.LeafScale;
    LocalPoint.z := LeafPolygonVertices[LeafShape, i, 2];

    with FTree.FModelParameters do
    begin
      { Scale the length of the leaf }
      LocalPoint.y := LocalPoint.y * LeafScale / Sqrt(LeafQuality);
      { Scale the width of the leaf }
      LocalPoint.x := LocalPoint.x * LeafScale * LeafScaleX / Sqrt(LeafQuality);
    end;

    GlobalPoint := AFrame.Local2Canonical(LocalPoint);
    Points.Add(GlobalPoint);
    LocalPoint.Destroy;

  end;
  FLeafs.Add(Points);
end;

procedure TStem.CreateFruit(const AFrame: TCoordinateFrame);
var
  i: Integer; // Loop variable
  Points: TList; // List of points
  LocalPoint: TPoint3d; // LocalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // local coordinates)
  GlobalPoint: TPoint3d; // GlobalPoint is used to hold the points
  // that are being created and added to
  // the points list of the section (in
  // global coordinates)
begin
  Points := TList.Create;
  { Depending on range use shape defined
    or simplest..    to less shapes+bigger leaves }
  { THIS is WHERE the Leaf Patterns would go
    Thus be able to make Shaped Globs of Leaves...
    But as Multiple .. Counts WRONG }
  for i := 0 to FruitPolygonVerticesNumber do
  begin
    LocalPoint := TPoint3d.Create;
    LocalPoint.x := FruitPolygonVertices[i, 0];
    { Get out of Branch Axis... invisible Stalk }
    LocalPoint.y := FruitPolygonVertices[i, 1] +
      FTree.FModelParameters.FruitScale;
    LocalPoint.z := FruitPolygonVertices[i, 2];

    with FTree.FModelParameters do
    begin
      { Scale the length of the leaf }
      LocalPoint.y := LocalPoint.y * FruitScale / Sqrt(LeafQuality);
      { Scale the width of the leaf }
      LocalPoint.x := LocalPoint.x * FruitScale * FruitScaleX /
        Sqrt(LeafQuality);
    end;

    GlobalPoint := AFrame.Local2Canonical(LocalPoint);
    Points.Add(GlobalPoint);
    LocalPoint.Destroy;
  end;
  FFruit.Add(Points);
end;

{ DOESN'T FUNCTION PROPERLY }
function TStem.AdjustLeafFrame(const ALeafFrame: TCoordinateFrame)
  : TCoordinateFrame;
{ Pre: ALeafFrame <> nil
  Ret: A leaf frame so that the orientation of the leaf is adjusted to
  optimize the available direct sun light and scattered sky light. }
var
  ThetaPosition: real;
  ThetaBend: real;
  PhiBend: real;
  Orientation: real;
  RotationMatrix: TGLMatrix;
  LocalX, LocalY: TVector3d;
  GlobalX, GlobalY: TVector3d;
  TempONB, TempONB2: TONB;
  TransformZ: TVector3d; // Vectors holding the unit y and
  // z vectors passed through the
  // current frame transformation
begin
  Result := TCoordinateFrame.Create;

  with ALeafFrame.Origin do
    ThetaPosition := ArcTan2(y, x);
  with ALeafFrame.ONB.w do
    ThetaBend := ThetaPosition - ArcTan2(y, x);

  RotationMatrix := TGLMatrix.Create;

  RotationMatrix.SetRotationMatrix(1 * ThetaBend, aZ);

  LocalX := RotationMatrix.Multiply(UnitX);
  LocalY := RotationMatrix.Multiply(UnitY);

  GlobalX := ALeafFrame.ONB.LCTransformMatrix.Multiply(LocalX);
  GlobalY := ALeafFrame.ONB.LCTransformMatrix.Multiply(LocalY);

  TempONB := TONB.CreateFromUV(GlobalX, GlobalY);

  with ALeafFrame.ONB.w do
    PhiBend := ArcTan2(Sqrt(x * x + y * y), z);

  LocalX.Destroy;
  LocalY.Destroy;
  GlobalX.Destroy;
  GlobalY.Destroy;

  TransformZ := TempONB.LCTransformMatrix.Multiply(UnitZ);
  Orientation := ArcCos(TransformZ.y);

  RotationMatrix.SetRotationMatrix(-Orientation, aZ);

  LocalX := RotationMatrix.Multiply(UnitX);
  LocalY := RotationMatrix.Multiply(UnitY);

  GlobalX := TempONB.LCTransformMatrix.Multiply(LocalX);
  GlobalY := TempONB.LCTransformMatrix.Multiply(LocalY);

  TempONB2 := TONB.CreateFromUV(GlobalX, GlobalY);

  LocalX.Destroy;
  LocalY.Destroy;
  GlobalX.Destroy;
  GlobalY.Destroy;
  TempONB.Destroy;

  RotationMatrix.SetRotationMatrix(PhiBend, aX);

  LocalX := RotationMatrix.Multiply(UnitX);
  LocalY := RotationMatrix.Multiply(UnitY);

  GlobalX := TempONB2.LCTransformMatrix.Multiply(LocalX);
  GlobalY := TempONB2.LCTransformMatrix.Multiply(LocalY);

  TempONB := TONB.CreateFromUV(GlobalX, GlobalY);

  RotationMatrix.SetRotationMatrix(Orientation, aZ);

  LocalX := RotationMatrix.Multiply(UnitX);
  LocalY := RotationMatrix.Multiply(UnitY);

  GlobalX := TempONB.LCTransformMatrix.Multiply(LocalX);
  GlobalY := TempONB.LCTransformMatrix.Multiply(LocalY);

  TempONB2 := TONB.CreateFromUV(GlobalX, GlobalY);

  Result.ONB.Assign(TempONB2);
  Result.Origin.Assign(ALeafFrame.Origin);
  Result.SetTransformationLists;

  { Clean up dynamically created objects }

  LocalX.Destroy;
  LocalY.Destroy;
  GlobalX.Destroy;
  GlobalY.Destroy;
  RotationMatrix.Destroy;
  TransformZ.Destroy;
  TempONB.Destroy;
  TempONB2.Destroy;
end;

initialization

{ oval, thin, triangle, (?3-lobe oak?), 3-lobe maple,
  5-lobe maple, 5-lobe oak, 3 leaflets, User Defined }
{ 7:0..6(types)  ?(points per leaf)  3:0..2  (x,y,z) }
{ LeafPolygonVertices: array of array of array of real; }
{ array[0..0, 0..8, 0..2] of real = }
SetLength(TreePolygonVertices, 6);
SetLength(TreePolyImageVertices, 6);
SetLength(TreePolygonVerticesNumber, 6);
TreePolygonVerticesNumber[0] := -1;
TreePolygonVerticesNumber[1] := -1;
TreePolygonVerticesNumber[2] := -1;
TreePolygonVerticesNumber[3] := -1;
TreePolygonVerticesNumber[4] := -1;
TreePolygonVerticesNumber[5] := -1;

SetLength(LeafPolygonVertices, 2);
SetLength(TreePolyImageVertices, 2);
SetLength(TreePolygonVerticesNumber, 2);
SetLength(LeafPolygonVertices[0], 9, 3);
{ (0, 0, 0) }
LeafPolygonVertices[0, 0, 0] := 0;
LeafPolygonVertices[0, 0, 1] := 0;
LeafPolygonVertices[0, 0, 2] := 0;
{ (-0.2, 0.1, 0) }
LeafPolygonVertices[0, 1, 0] := -0.2;
LeafPolygonVertices[0, 1, 1] := 0.1;
LeafPolygonVertices[0, 1, 2] := 0;
{ (-0.3, 0.5, 0) }
LeafPolygonVertices[0, 2, 0] := -0.3;
LeafPolygonVertices[0, 2, 1] := 0.5;
LeafPolygonVertices[0, 2, 2] := 0;
{ (-0.2, 0.9, 0) }
LeafPolygonVertices[0, 3, 0] := -0.2;
LeafPolygonVertices[0, 3, 1] := 0.9;
LeafPolygonVertices[0, 3, 2] := 0;
{ (0, 1, 0) }
LeafPolygonVertices[0, 4, 0] := 0;
LeafPolygonVertices[0, 4, 1] := 1;
LeafPolygonVertices[0, 4, 2] := 0;
{ (0.2, 0.9, 0) }
LeafPolygonVertices[0, 5, 0] := 0.2;
LeafPolygonVertices[0, 5, 1] := 0.9;
LeafPolygonVertices[0, 5, 2] := 0;
{ (0.3, 0.5, 0) }
LeafPolygonVertices[0, 6, 0] := 0.3;
LeafPolygonVertices[0, 6, 1] := 0.5;
LeafPolygonVertices[0, 6, 2] := 0;
{ (0.2, 0.1, 0) }
LeafPolygonVertices[0, 7, 0] := 0.2;
LeafPolygonVertices[0, 7, 1] := 0.1;
LeafPolygonVertices[0, 7, 2] := 0;
{ (0, 0, 0) }
LeafPolygonVertices[0, 8, 0] := 0;
LeafPolygonVertices[0, 8, 1] := 0;
LeafPolygonVertices[0, 8, 2] := 0;

LeafPolygonVerticesNumber[0] := 8;

{ oval, thin, triangle, (?3-lobe oak?), 3-lobe maple,
  3-5-lobe maple, 5-lobe maple, 5-lobe oak, 3 leaflets,
  User Defined  This is just a start... a bad one }
SetLength(LeafPolygonVertices[1], 2, 3);
{ User Defined   (0, 0, 0) }
LeafPolygonVertices[1, 0, 0] := 0;
LeafPolygonVertices[1, 0, 1] := 0;
LeafPolygonVertices[1, 0, 2] := 0;
{ User Defined   (0, 0, 0) }
LeafPolygonVertices[1, 1, 0] := 0;
LeafPolygonVertices[1, 1, 1] := 0;
LeafPolygonVertices[1, 1, 2] := 0;
(* LeafPolygonVerticesNumber: array{[0..0]} of integer{ = (9)}; *)
LeafPolygonVerticesNumber[1] := -1;

finalization

SetLength(LeafPolygonVertices, 0);
SetLength(LeafPolyImageVertices, 0);

SetLength(FruitPolygonVertices, 0);
SetLength(FruitPolyImageVertices, 0);

SetLength(TreePolygonVertices, 0);
SetLength(TreePolyImageVertices, 0);
SetLength(TreePolygonVerticesNumber, 0);

end.
